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STATUS REPORT - JUNE 1974 


NASA AMES GRANT NGR 37-008-003* 

During the subject reporting period the research discussed herein was advanced by 
a graduating mechanical engineering senior, Daniel W, Drago, and the principal investi- 
gator. As a consequence of the previously unprogrammed investigation of some interesting 
aspects of radiative transfer calculations for non-unity index of refraction (discussed 
below), a request for extension of the completion date of the present research to August 31, 
1974 was requested and recently approved. Mr. Drago is remaining at the University as a 
graduate student and will be participating through the grant termination date. 

In the present reporting period two papers have been submitted for publication and 
another paper published. As reported in the last status report, a paper by R. S. Reddy and 
the principal investigator had been submitted to and accepted for publication by the AIAA 
Journal . This paper, which has since appeared in the Journal (reference 1), provides an 
approximate analytic solution for the unsteady radiative heating of a highly scattering 
medium under constant heat flux boundary conditions. 

A more recent paper prepared under this grant (reference 2) explores the accuracy 
of two-flux methods (such as that employed in reference 1) in representing the reflectance, 
radiative flux and radiative flux divergence of highly scattering materials. The influence 
of anisotropic scattering on these parameters was also considered. This work, not accepted 

for presentation at the 1974 Heat Transfer and Fluid Mechanics Institute, has been revised 
and submitted for publication in the AlAA Journal. In addition to defining conditions of 

*The NASA Technical Officer for this grant is Dr. Phillip R. Nachtsheim, NASA 
Ames Research Center, Moffett Field, CA 94035. 
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applicabilty of approximate methods, this work clarifies relations between transfer 
equation solutions and approximate methods, and presents new results for anisotropic 
scattering . 

Still more recently, while evaluating the capabilities of the NANITERAD 
program in dealing with radiative transfer In media with non-unity refractive index, 
solution intensities were found to vary erratically with quadrature order. Through an 
approximate discontinuity model and transfer equation solutions we have shown that, 
for a fixed quadrature, the quadrature accuracy (and hence transfer equation solution 
accuracy) varies erraficaHy with refractive index rather than smoothly as implied in the 
literature (reference 3 for Instance). The discontinuity model shows that under certain 
conditions the accuracy of Gaussian quadrature may decrease with increasing quadrature 
order rather than increase as found with continuous distributions. The observed error 
pattern is due to the relative position of the critical angle for total internal reflection 
and its neighboring quadrature directions. It was found that the errors induced by the 
discontinuity are amplified in an iterative solution due to the scattering term in the trans- 
fer equation. A published solution (reference 3) for refractive index of 1 .4 was shown to 
exhibit a large discontinuity error as predicted by the model . As a result of this study we 
are revising our calculation procedures to Incorporate separate quadratures within and 
outside the critical angle for total internal reflection. For a continuous distribution the 
accuracy of the numerical integration is reduced by using several quadrature formulae in 
place of a single quadrature application with the same total number of directions. How- 
ever for the non-unity index cases, this loss may be far exceeded by the gain in accuracy 
achieved by discontinuity error avoidance. A discussion of the errors Incurred by the use 
of single quadrature has been presented in a paper submitted to the ASME Journal of Heat 



Transfer (reference 4). We have been advised thaf the general reaction of the ASME 
reviewers was favorable. The paper has been resubmitted incorporating certain revisions 
as requested. 

The formation of a minimum in the magnitude of reflected intensities near the 
rear surface of an absorbing and scattering slab is considered in Appendix A of this 
report. The intensities decrease, stablize and increase as the front surface is approached 
This unexpected phenomenon was first considered to occur because of a programming 
error. A simple analytical distribution model applied to the transfer equation however 
verified the existence of a dip under certain conditions. Model and transfer equation 

results are also compared in the Appendix. 

A study of the most suitable quadrature formulae combination with a limited 
number of quadrature directions for non-unity index of refraction has been completed 
and appears as Appendix B to this status report . This study shows that the most accurate 
method of formulae application for cases with significant discontinuities uses three sepa- 
rate quadratures over the range one each within the positive and negative 

regions of total Internal reflection and the third between the other two. The appendix 
also provides tabulated error results for various single and multiple formulae applications 

The NANITERAD program is at this time regarded as fully operational . A 
number of cases have shown consistent and reasonable behavior for non-unity refractive 
index although the results available in the literature for comparison are quite limited. 
The most convincing check was with a fifth order quadrature calculation given in 
reference 4. The detailed radiation field agreed with the NANITERAD calculations 
within reading errors dictated by the small size of the figure of the reference. 
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APPENDIX A 

THE CONDITIONS FOR INTENSITY EXTREMA 
While exploring various aspects of the transfer equation results, the radiative 
intensities reflected from the rear surface were thought to monotonicaHy increase or 
decrease towards the front surface. Occasionally the reflected intensities were found 
to decrease near the rear boundary, stabilize, and increase as the front boundary is 
approached, forming a distinctive dip. A simple model of the transmitted and reflected 
intensifies is applied to the transfer equation, identifying conditions which allow the 
minimum intensity to exist between the boundaries. 

The intensities at any distance y from the front surface are modeled as: 

f I -(y), n<0 

(Al) 

+ (y),H>o 


+/ 


where 1 (y ) Is the transmitted intensity and I (y ) the reflected intensity. 

By using (Al) in the transfer equation while neglecting the emission term, and 
considering the variation of 1 (y ) m the direction corresponding to -| li*| the slope of 
I ”(y ) becomes (for 1 H*| +0): 


d i ~(y) 

dy 


= [qj:> /K + sj-il 

W\ |Pfy)' { ? 2_). 


As a special case, with pure absorption (A2) reduces to: 
d I ~(y) = JKJ“(y) 

dy I n*| 


(A2) 


(A3) 


In this case the minimum occurs at the front surface since the slope is always positive or 


zero . 
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In a second special case, with pure scattering (A2) yields: 

'(V) 


di -(/) = i (y) / s \ 

dy |u*| ' 2 
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Since the right hand side of A4) is always negative or zero for a flux incident at y - 0, 
the minimum of 1 “(y) is at the rear boundary. Thus in either the pure absorption or pun 

scattering case the dT pf.- may not occur. 

For non-vanishing scattering and absorption coefficients an extremum of I (y) 
occurs when the right side of (A2) vanishes at some point within the medium. At that point 
(A2) indicates a condition for an extremum to be: 


1 ~(y ) = 1 

TW) 


2 | + i . 


(A5) 


It may be shown using (A2) and a similar equation for the slope of l + (y) that the 
second derivative of I “(if) is always positive or zero; 


,,d 2 r,(y) d f(y) s. «*! *(-r> ,S % 
■pf = (K + T ) “ <2 } 


dy : 


dy 


dy 


= [ l “(y) (K + 2 ) “ I + (y) (2 )1 + 2) 

- l-i + (y) (k + 7) + 1 “(y) (f) ] (f) 


or 


d 2 I ~(y) = l~(y)(K 2 +KS) 

d y2 1 pi*| ^ 


(A6) 


Since the right side of (A6) is always non-negative, this precludes the possibility of an 
extremum being a maximum rather than a minimum. Therefore, if the slope of 1 ~(y ) at 
the rear boundary (/ - y ‘ ) is shown to be negative or zero, a dip will not occur. This 
slope at the rear surface may be calculated to determine the conditions necessary for the 
slope to be positive, and thus to provide a third condition for the formation of the dip. 
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The intensities at the rear surface are modeled as: 


Cr l + ( y ) , u<0 

I W = + ° 

( i (y o ) , H>0 

where R„ is the' rear surface reflectivity. 

D 

Evaluating (A2) at y Q using (A7) produces the slope of I (y Q ): 


d I~(y 0 ) = i (y p ) 

dy |M.*| 




(A 7) 


(A8) 


A positive slope of 1 "(Vq) exists when: 
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Thus for the mode! given by (Al) and A7), three conditions must be satisfied for a 
dip to exist: non-zero scattering and absorption coefficients; a sufficiently large R^ 
consistent with (A9); and a small 1 “(y) consistent with (A5). The latter condition is 
satisfied only when the optical thickness is large enough for l (y) to obtain its minimum 
within the boundaries. 

The criterion fora minimum from (A9) is evaluated using transfer equation solutions 
in Figure A 1 . This figure shows the conditions under which minima were found in the solu- 
tions. Due to the deviation of the actual intensify field at the rear surface from that 
modeled in (Al), the condition for a dip based on (A9) has been compared to transfer 
equation solutions demonstrating a dip in a majority of the quadrature directions. The 
data from the solutions has been divided info three categories: 1) None of the eight 
Gaussian quadrature directions contain a dip; 2) At least one direction displays the mini- 
mum; and 3) Five or more directions have a dip. As seen in Figure Al, the minimum Rg 
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necessary to produce dips in five directions closely follows the trend of the pattern set by 
condition (A9). 

In Figure A2 the optical depth, T, of the minimum from the rear surface is clearly 

seen to be a function of H . A similar pattern of minima is evident for each case studied. 

This pattern may be qualitatively explained by referring to equations (A2) and (A6). 

Assuming a directionally independent value for the reflected beams at the rear boundary 

and an actual R g larger than the R B _ cr | t | ca j from (A 9), equation (A2) indicates the 

reflected beams diminish most rapidly parallel to the boundary. Equation (A6) is even 

2 

more dependent on M- than (A2) due to its inverse relationship to . As such, a plot of 
the intensity of a beam directed close to parallel to the boundary would demonstrate much 
larger curvature, and hence would stabilize and begin to Increase within substantially 
less distance from the rear surface, than a beam more nearly normal to the boundary. 

A second observation from Figure A2 is that the depth of the minimum from the 
rear surface in each Gaussian direction generally increases as the ratio of the actual Rg 
to the model Rg - cr Ttr ca l * ncreases / a P°^ ern further explored in Figure A3. Here the 
depth of the minimum in the direction closest to normal to the boundaries of the medium 
is plotted against the ratio of the Rg’s. The distance of the minimum from the rear surface, 

T , rapidly increases as the rear surface reflectivity is increased. Lines with constant 
values of k/s have generally increasing slope as the value of k/s decreases, with the 
smallest plotted value of k/s producing a nearly vertical line. The actual intensity field 
becomes more closely approximated by the model in (AT) as k/s approaches zero, and as 
such increases the accuracy of prediction based on (A9). Thus, this study shows that minima 
should be expected to occur over a wide range of conditions for highly scattering media. 
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APPENDIX B 


THE ACCURACY OF GAUSSIAN QUADRATURE COMBINATIONS 

As discussed in [4], a significant error may be generated in the numerical approxi- 
mation of the scattering integral in the transfer equation when a single quadrature is used 
for media with refractive indices greater than one. The error arises from the presence of 
a discontinuity in the intensity field at the critical angle for total internal reflection due 
to a Fresnel boundary. Figure B1 shows the computed values of transmission across such a 
boundary when the quadrature is applied over three different regions. The first and second 
applications, which are indifferent to the presence of the discontinuity, result in large and 
erratic errors for small to medium quadrature orders. The third application restricted to the 
region of transmitted flux only, demonstrates high accuracy even for low quadrature order 
and smoothly converges to the exact value. 

The overall accuracy in the transfer equation solutions when using single quadrature 
applications has been shown to follow a similar pattern to that found with the first two 
applications used in Figure B1 [4]. In this appendix the question of which combination of 
Gaussian quadratures is most suitable for radiative transfer calculations Involving discon- 
tinuities at the critical angle is considered, with the overall medium reflectance providing 
a measure of combination quality. The quadrature combinations are applied such that the 
quadratures selected are bounded by the critical angle and have a combined total of no more 
than sixteen quadrature points. Three regions ofM-were defined for the application of 
.quadrature formulae to the scattering integral: 
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Region 1* -1 <^i<ii c 

Region II: -\l c < |-l< M c (81) 

Region ill : M-c ^ ^ — 1 

where \i c is the cosine of the critical angle of total internal reflection: 

H C =y 1 - (l) 2 (B2) 

and N the refractive index of the medium. Typically, identical quadrature orders were 
applied to Regions 1 and III, with Region 11 utilizing the remainder of the sixteen quadrature 
directions. 

The notation used here for the quadrature formulae application Is: 

X/VA < B3 ) 

where X is the quadrature order of Regions I and 111, and Y is the order of Region 11 . For 
example, fourth order quadratures may be used in Regions I and 111 with eighth in Region II, 
as expressed by 4/8/4. 

Reflectance results for various quadrature combinations used in approximating the 
transfer equation scattering integral are presented in Table Bl. Selected combinations were 
run for two refractive Indices to insure that the relative magnitudes of the reflectance values 
were not significantly affected by a refractive index change. The same reflectance pattern 
did appear in both cases, with S/3/3/5 (Region II was divided at M = 0 with two quadratures 
applied) and 4/8/4 nearly identical, 6/4/6 lower, and 5/6/5 between. As predicted In [4], 
tenth quadrature was found to have a large error at N = 1.4 while sixteenth was relatively 
accurate. 

As shown in Figures B2 and B3, a reference reflectivity was determined from three 
sequenced test values found for each refractive index case. These values came from 
quadrature combinations of 2/4/2, 3/6/3, and 4/8/4, with the values of reflectivity 
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plotted as a function of the inverse number of total quadrature directions. The reference 
value was Found by extrapolating a line faired on the three test points to the zero value of 
the abscissa. 

Study of these figures and the Table show that both 5/6/5 and 5/3/3/5 are within 
.12% of the reference values at both refractive indices, with the average percentage 
error smallest for the 5/6/5 combination. Thus the 5/6/5 application was ultimately 
selected for use in the transfer equation solutions for cases with non-unity index of 
refraction. A single quadrature application of sixteenth order was retained for the N — 1.0 


cases . 




TABUS B1 

TABLE OF REFLECTANCE VALUES AND ERRORS 
FOR SEVERAL QUADRATURE COMBINATIONS 



OJ = 1.0 

0* K 




'r a = i.o 

I - 0 




Bb«.Q3 
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QUADRATURE 

N = 

1.4 

N« 

1.2 

-COMBINATION 

R 

E$) 

H 

E(3) 

4/8/4 

.290044 

•.0807 

.344580 

-.1570 

5/3/315 

*290001 

-.0659 

.344429 

-.1131 

5/6/5 

*289466 

.1187 

*344042 

-.0006 

6/4/6 

*288170 

.5659 

.343205 

.2427 

10 

*220740 

23.83 
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*288401 

.4862 


— 

- 2/4/2 

*290637 

-.2854 

.345852 

-.5267 

3/6/3 

*290197 

-.1335 

.344879 

-.2439 

Reference 

*289810 

0 

-*344040 

0 


Value 
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